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SUMMARY 

The operation of the INS3D code, which computes steady-state solutions to the incom- 
pressible Navier-Stokes equations, is described. The flow solver utilizes a pseudocompressibility 
approach combined with an approximate factorization scheme. This publication describes key 
operating features to orient new users. This includes the organization of the code, description 
of the input parameters, description of each subroutine, and sample problems. Details for more 
extended operations, including possible code modifications, are given in the Appendix. 


INTRODUCTION 

This publication is intended as a reference guide for users of the INS 3D code (ref. 1). 
Many recent advances in computational fluid dynamics (CFD) have been made in conjunction 
with compressible flow computations. Therefore, the code INS3D is structured to take advantage 
of fast algorithms developed for compressible flows. The primary features can be summarized as 
follows. 

(1) A pseudocompressibility approach is used in which a time derivative of pressure is added 
to the continuity equation, which together with the momentum equations form a set of 
four equations with pressure and velocity as the dependant variables (refs. 2, 3, 4). 

(2) The coordinates are transformed for general three-dimensional applications. 

(3) This set of equations is solved using an approximate-factorization scheme (refs. 5, 6). 

(4) Computational efficiency is achieved using a diagonal algorithm (ref. 7). A block tri di- 
agonal option is also available. 

(5) When a steady-state solution is reached, the modified continuity equation will satisfy 
the divergence-free velocity field condition. 
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The code has been verified by computing fundamental fluid dynamics problems such as 
the flow through a channel (ref. 1), flow over a backward-facing step (ref. 8), and flow over a 
circular cylinder (refs. 1 and 8). Three-dimensional cases include flow over an ogive cylinder(ref. 
9), flow through a rectangular duct(ref. 1), wind-tunnel-inlet flow (ref. 10), cylinder- wall juncture 
flow and flow through multiple posts mounted between two plates (refs. 11 and 12). The most 
striking application of the present code is the flow simulation in the power head portion of the 
Space Shuttle Main Engine (refs. 13 and 14). 


GENERAL USER NOTES 

This section contains some important features of INS3D which will be useful for the 
first-time user. The arrays in the program are dimensioned using a single index I, which is a 
function of the three single-direction indices. 

I = J -h (K-1)*JMAX + (L-1)*JMAX*KMAX 

Here, J is the index in the ^-direction, K is the index in the » 7 -direction, and L is the index in 
the ^-direction. The grid point locations are stored in the single index arrays X, Y, and Z. The 
primitive variables are stored in the array Q, where 

Q(1,I) = pressure 
Q(2, 1) = u velocity component 
Q(3,I) = V velocity component 
Q(4, 1) = w velocity component 

To run a specific geometry, the user will have to change a number of different subroutines. 
These subroutines are GRID, IC, BC, and SMOOTH. The details of each of these changes are 
included here. 

The grid coordinates for the problem should be loaded into the X, Y, and Z arrays in 
subroutine GRID. This subroutine is called once at the beginning of the program execution. If 
the grid can be generated algebraiccilly, this can be done in this subroutine. If the coordinates 
are generated with some other program, then the coordinates merely have to be read in and 
stored in the X, Y, and Z arrays. In defining the geometry for the problem, the coordinates 
should be normalized by the characteristic length of the problem because the program assumes a 
characteristic length of unity when it uses the input Reynolds number. 

The subroutine IC should be changed to set the initial conditions for the the specific 
problem. It is usually adequate to set the pressure everywhere to unity, and to set the velocity 
components to freestream conditions everywhere except on no-slip walls, where zero velocity 
should be specified. The velocity field should also be normalized with the characteristic velocity 
for the problem because the program assumes a characteristic velocity of unity when using the 
Reynolds number. 

The third subroutine which needs to be changed is BC, where the explicit boundary 
conditions are applied. By default, the code will update all the interior points during one iteration. 
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and then call subroutine BC. See the Numerical Algorithm section in the appendix for more details 
on the proper application of boundary conditions. 

Changes may also have to be made to the subroutine SMOOTH, which adds the artificial 
dissipation to the right side of the equations. In general, fourth-order explicit smoothing is added 
at all the interior grid points. Stability of the calculations is enhanced if second-order smoothing 
is added at the interior points adjacent to the inflow and outflow boundaries. Sometimes when 
stability problems occur near severe geometries or grid areas, the use of second-order explicit 
smoothing will help to smooth out the disturbance. 

Another change of primary importance is to correctly set the dimensions of the arrays 
for the specific problem. This is easily done with the use of one parameter statement which is 
included in almost all the subroutines. Two parameters are included in this statement, JKLMAX 
and IMX. The first parameter JKLMAX should be set to the majdmum of JMAX, KMAX, or 
LMAX. The second parameter IMX should be set equal to the total number of grid points, that is, 
the product of JMAX, KMAX, and LMAX. All of the arrays are passed between the subroutines 
with the use of common blocks. 

Periodic cases can be handled if the problem is set up so that the periodic direction 
corresponds with the ?;-direction. The K=KMAX and K=1 planes must be adjacent to each 
other and are treated as interior points in this case. 

Two-dimensional problems can be handled by setting JMAX, KMAX, or LMAX to three. 
For instance, if the two-dimensional problem is chosen to lie in the JK plane in computational 
space, then LMAX should be set to three. Assuming that the physical space is set up to lie in the 
xy plane, the grid should be duplicated in three z = constant planes. The dependant variables 
should be initialized to have identical values in the L=l, 2, and 3 planes. No calculations will 
be done in the direction which has a maximum dimension of three, and calculations will only be 
done on the middle of the three planes. A loop should be put in subroutine BC which copies the 
data from the middle plane to the two outer planes after each iteration. 
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VARIABLES IN NAMELIST DATAIN 


The input for the program is read in with the use of a namelist. Each of the variables 
in this namelist are described here. This is done in tabular form; the entries are the name of 
the variable, the possible range of acceptable values, itsi default value, and a description of the 
variable. 


Variable 

Range 

Default 

Description 

BETA 

0.1 to 50.0 

5.0 

This is the artificial compressibility parameter. The 
artificial speed of sound is a function of this param- 
eter, and the minimum value this should take is de- 
pendent on the geometry. The majdmum value is 
bounded because of the error introduced by the ap- 
proximate factorization. For more details, see the 
pseudocompressibility section in the Appendix. 

DISKOUT 

True or False 

False 

This is a logical variable which controls whether sub- 
routine DISKIO is called at the end of the program 
to write the Q array to disk. The array is written 
out to unit 51. 

DTAU 

0.0001 to 0.1 

0.05 

This is the time step used at each iteration. If a 
stability problem occurs, this parameter should be 
reduced. 

DXDT 

— oo to oo 

0.0 

This is a parameter which is used to specify transla- 
tion of the grid in the x-direction. It is used when 
the time varying metrics are computed. 

ENDACC 

0 or 1 

1 

This is an integer which controls the spatial differ- 
encing of the metrics at points next to nonperiodic 
boundaries. If ENDACC=1, first-order-accurate one- 
sided differencing is used. Otherwise, second-order- 
accurate one-sided differencing is used. 

IBLKDIA 

1 or 2 

2 

1 

1 

This integer determines whether the code will use a 
standard block algorithm for the inversion of the ap- 
proximately factored matrices or a diagonalized al- 
gorithm. For IBLKDIA = 1, the standard algorithm 
is used; for IBLKDIA = 2, the diagonal algorithm is 
used. The standard algorithm inverts 4-by-4 block 
tridiagonal matrices whereas the diagonal algorithm 
inverts scalar tri diagonal matrices. The diagonal al- 
gorithm is about three times faster for a single time 
step than the standard algorithm, and has the same 
convergence characteristics. 
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Variable 

Range 

Default 

Description 

IMPSMO 

2 or 4 

2 

This integer determines the order of the implicit arti- 
ficial dissipation when the diagonal algorithm is used. 
When set equal to 2, second-order dissipation is used, 
when set equal to 4, fourth-order dissipation is used. 
Not much difference has been seen between the two 
different orders. The fourth order is slightly more 
expensive as it requires the inversion of scalar penta- 
diagonal matricies. 

lORTHO 

1 or 2 

2 

This integer affects the assumptations in computing 
the viscous terms. For I0RTH0=1, the nonorthogo- 
nal contributions to the viscous terms are included in 
the right-hand side of the equations by calling subrou- 
tine VISRHS2. For I0RTH0=2, the viscous terms 
are simplified by assuming an orthogonal grid. 

IPRNT 

1 to NTMAX 

1 

This integer controls the frequency with which the 
convergence history is written. Convergence data is 
written every IPRNT*^ iteration. 

ISTART 

0 or 1 

0 

This variable is an integer which controls the restart 
capability of the code. If ISTART = 0, then the 
code starts from the conditions specified in subrou- 
tine IC. If ISTART = 1, then the code calls subrou- 
tine DISKIO which attempts to read in NT, JMAX, 
KMAX, LMAX and the Q array from unit 20. The 
grid is always computed in subroutine GRID regard- 
less of the value of ISTART. 

JMAX 

> 3 

52 

This integer is the maximum number of points used 
in the ^-direction. 

KMAX 

> 3 

63 

This integer is the maximum number of points used 
in the ?;-direction. 

KPERI 

0 or 1 

1 

This integer allows one to specify a periodic boundary 
between the K=KMAX and K=1 planes. If KPERI 
= 1, the boundary is assumed to be periodic and the 
points at K=1 and K=KMAX are treated as interior 
points. If KPERI = 0, these points are treated as 
normal boundary points. 

LMAX 

> 3 

35 

This integer is the maximum number of points used 
in the ^-direction. 
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Variable 

Range 

Default 

Description 

NPRNT 

> 1 

50 

This integer controls the frequency with which flow 
solution output is to be written during the iteration 
process. Every NPRNT**' iteration subroutine OUT- 
PUT is called, which prints out the values for Q, X, 
Y, Z, Cp, and the Jacobian for specified indices. 

NTMAX 

> 1 

100 

This integer is the number of iterations the program 
will run. 

PL0T3D 

True or False 

False 

This logical variable determines whether unformat- 
ted files to be used with the plotting program known 
as PLOT3D will be written at the end of the com- 
putation. If the variable is TRUE, the grid file is 
written to unit 30, and the Q file is written out to 
unit 31. 

REYNUM 

0.0 

1000 

This is the Reynolds number based on unit charac- 
teristic length and unit characteristic velocity. 

RLl 

0.0 or 1.0 

1.0 

This variable is used for a geometry in which the 
volume of the grid cells becomes zero at L = 1. This 
can occur with an 0-grid whose innermost circle has 
a radius of zero. This variable then represents the 
radius at the L = 1 circle, and is used only when it 
is set equal to zero. 

SMU 

0.0 

0.1 

This real number is the smoothing coefficient used 
for the explicit dissipation added to the right-hand 
side. This should be one-half to one-third the value 
of SMUIM. 

SMUIM 

0.0 

0.3 

This real number is the smoothing coefficient used for 
the implicit dissipation added to the left-hand side. 
This should be two to three times the value of SMU. 

SMUPRS 

0.0 

1.0 

This real number is a coefficient which multiplies the 
smoothing terms added to the right-hand side of the 
pressure equation in addition to SMU. It may become 
necessary to reduce this coefficient as the solution 
converges so that the divergence of the velocity will 
continue to converge toward zero. For more details, 
see the section on the effects of smoothing on pressure 
in the pseudocompressibility section of the Appendix. 
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DESCRIPTION OF ROUTINES AND CODE STRUCTURE 


In this section, a brief description of the subroutines found in INS3D is given. A flow 
chart of the main code and of subroutine STEP is then presented. 

Program MAIN 

Program MAIN handles the main control of the program. It flrst calls routine INITIA 
and then enters the main time-iteration loop. Inside this loop, the routines STEP, BC, and VISCT 
are called. If printing of the solution is requested for the current iteration, then OUTPUT is also 
called. After the iteration process is finished, MAIN also handles the writing of the solution to 
disk by calling DISKIO and P3D0UT if requested. 

Subroutine AMATRX 

Subroutine AMATRX computes the Jacobian matrix Ai and loads it into the D array in 
common block DUDD. It is computed and stored for all the points along a single grid line (two of 
the three j,k,l indicies held constant). The proper metrics must be loaded into the metric storage 
array depending on whether A,R,or C is being computed. The Jacobian matrix is used to All 
the implicit side matrices. 


Subroutine BC 

Subroutine BC enforces the boundary conditions explicitly. It is called at the end of 
each iteration, after all the interior points have been updated to the new time-step. This routine 
must be changed for each specific geometry run. 

Subroutine BTRI4 and PBTRI4 

Subroutines BTRI4 and PBTRI4 invert a system of block tridiagonal matrices where 
the block size is 4 by 4. The system of equations must be loaded into the matrices in the BTRI 
common block. The solution is overwritten onto the forcing-function array F. BTRI4 handles the 
nonperiodic case and PBTRI4 handles the periodic case. 

Subroutine COMET 

Subroutine COMET computes and stores the elements of the first two eigenvectors of 
the Jacobian matrix. These quantities are used in computing the similarity transformation which 
diagonalizes the Jacobian matrix. This routine must be called prior to calling TK or TKINV. 

Subroutine DISKIO 

Subroutine DISKIO is used to read or write the flow field solution (the array) from the 
disk. Data is read from fortran unit 20, and is written to fortran unit 51. 
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Subroutine ETAINV, XIINV, and ZETAINV 


Subroutines ETAINV, XIINV, and ZETAINV load the implicit side matrices for each 
of the three different sweeps of the diagonal algorithm. The eigenvalues of the Jacobian matrix 
are computed and differenced to form the elements of the diagonal matrices. The implicit viscous 
terms are added and then either second- or fourth-order artificial dissipation terms are added. The 
system is then inverted by the appropriate scalar inversion routines. Since the first two eigenvalues 
are identical, the left-hand sides of the first two equations are identical; only the forcing function 
is different. Therefore these two are inverted together, and the last two equations are inverted 
separately. 


Subroutine FILTRJ, FILTRK, FILTRL, AND PFILTK 

Subroutines FILTRJ, FILTRK, FILTRL, and PFILTK load the implicit-side matrices 
for each of the three different sweeps of the block algorithm. PFILTK is used for the periodic 
case. These routines load the matrices with the difference of the Jacobian matrices and add the 
implicit viscous terms as well as second-order artificial dissipation terms. The system is then 
inverted using BTRI4 or PBTRI4. 


Subroutine GRID 

Subroutine GRID is used to compute the computational grid and load the X, Y, and Z 
arrays. GRID is called at the beginning of every run whether or not it is a restart. This routine 
must be changed for each specific geometry run. 

Subroutine IC 

Subroutine IC is called at the beginning of all nonrestart runs and is used to initialize 
the flow-field- array Q and the eddy-viscosity-array VNUT. The desired starting conditions must 
be coded for each specific geometry run. 

Subroutine INITIA 

Subroutine INITIA initializes some variables and arrays and reads in the input namelist 
to obtain the user-specified parameters. It calls routines GRID, JACOB, and then calls DISKIO 
if the run is a restart or IC if it is not. 

Subroutine JACOB 

Subroutine JACOB computes the Jacobian of the transformation for the given grid and 
stores it in array DJ. 
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Subroutine LUDEC and LUSOL 


Subroutines LUDEC and LUSOL perform L-U decompostions and solutions for a 4-by-4 
matrix. These subroutines are called by subroutines BTRI4 and PBTRI4. 

Subroutine METRIC 

Subroutine METRIC computes the metrics for a constant grid line and stores them in 
the metric arrays in common block METl. 

Subroutine OUTPUT 

Subroutine OUTPUT prints out the flow variables, the x, y, and z location, the Jacobian 
and the pressure coefficient for points along a given line. 

Subroutine P3DOUT 

Subroutine P3D0UT writes out the flow solution and the grid point variables in a format 
compatible with the PL0T3D plotting program. The grid-file is written to Fortran unit 30 and 
the q-flle is written out to Fortran unit 31. 

Subroutine PENTA, PENTA2, PENTAP, and PENTAP2 

Subroutines PENTA, PENTA2, PENTAP, and PENTAP2 invert a scalar pentadiagonal 
system of equations. A whole plane of equations are inverted at once and must be loaded into the 
arrays in the PNTD common block. A whole plane of data is done at once so that the loops can 
be fully vectorized in the nonsweep direction. These routines are called by the diagonal algorithm 
routines when fourth-order implicit artificial dissipation is used. 

Subroutine RHS 

Subroutine RHS computes the right-hand side of the governing equations by computing 
the spatial difference of the flux vectors at each point. Artificial dissipation is added by calling 
subroutine SMOOTH and the right-hand side viscous terms are added by calling subroutine 
VISRHS and subroutine VISRHS2. The right-hand side of the equations is stored in array S. 

Subroutine SMOOTH 

Subroutine SMOOTH adds the explicit fourth-order artificial dissipation terms to the 
right-hand side in the S array for all interior points. At most of the points next to the boundaries, 
a shifted fourth-order stencil is used. At points next to the inflow and outflow boundaries, second- 
order dissipation should be used to help maintain stability. This routine should be changed for 
each new geometry run. 
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Subroutine STEP 


Subroutine STEP marches the solution forward in time one iteration. It first calls RHS, 
and then successively calls three implicit-sweep routines. Depending on the value of IBLKDIA, 
either the block algorithm routines are called or the diagonal algorithm routines are used. After 
the inversion process is finished for all three sweeps, the fiow variables at all interior points are 
updated to the new time-step. 


Subroutine TK and TKINV 

Subroutines TK and TKINV load the eigenvector matrix (TK) and its inverse (TKINV) 
and multiplies them by the forcing function stored in the S array. This is done for a single grid 
line at a time. The proper metrics must be loaded in the common block MBUF , and the elements 
of the first two eigenvectors must be loaded by the subroutine COMET. These routines are used 
in the diagonal inversion process. 

Subroutine TRI, TRI2, TRIP, and TRIP2 

Subroutines TRI, TRI2, TRIP, and TRIP2 invert a scalar pentadiagonal system of 
equations. A whole plane of equations are inverted at once and must be loaded into the arrays 
in the PNTD common block. A whole plane of data is done at once so that the loops can be 
fully vectorized in the nonsweep direction. These routines are called by the diagonal algorithm 
routines when second-order implicit artificial dissipation is used. 

Subroutine VISRHS and VISRHS2 

Subroutines VISRHS and VISRHS2 compute the right-hand-side viscous terms and adds 
them to the S array for all the interior points. The viscous terms for an orthogonal mesh are 
computed in VISRHS, and the nonorthogonal terms are computed in VISRHS2. 
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INS3D FLOWCHART 



( STOP ) 
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SUBROUTINE STEP FLOWCHART 
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SAMPLE CALCULATIONS 


In this section three sample calculations are presented. The first example is a two- 
dimensional calculation of flow over a backward-facing step, and the second is the flow around a 
two-dimensional circular cylinder. The third case is flow around a three-dimensional post mounted 
normal and between two flat plates. Both the diagonal algorithm (IBLKDIA = 2) and the block 
algorithm (IBLKDIA = 1) were used for the backward-facing step, and the diagonal algorithm 
was used for the last two cases. The calculations were carried out on a Cray XMP-48. The 
average computer time on this machine used by the code for each case is approximately 3 x 10~® 
secs per grid point per iteration for the diagonal algorithm. The block algorithm is between 2.5 
and 3 times slower than this for most cases. 

Backward-Facing Step 

The flow over a backward-facing step with a two-to-one expansion ratio was computed 
for a laminar range of Reynolds numbers (based on twice the step height) between 100 and 800. 
A 65-by-33 algebraically-generated H-grid was used with the inflow boundary placed at the step. 
The grid was clustered near all the no-slip walls using an exponential stretching function. 

The relevant parameters used in this calculation are given in the following list. 


BETA 

= 5 

KPERI 

= 0 

DTAU 

= 0.1 

LMAX 

= 65 

ENDACC 

= 1 

SMU 

= 0.1 

lORTHO 

= 2 

SMUIM 

= 0.3 

JMAX 

KMAX 

= 3 
= 33 

SMUPRS 

= 0.1 


The length of the separation bubble behind the step was calculated for each Reynolds 
number and these results are compared to the experimental results of Armaly et al. (ref. 15). 
This comparison is shown in Figure 1, which shows good agreement between the experiment 
and the calculations for both the diagonal zilgorithm and the block algorithm. The convergence 
history is shown in Figure 2. Here, the rms of the change in the flow variables (RMSDQ) at each 
time-step is plotted versus iteration number. Three cases are plotted; (1) the block algorithm, 
(2) the diagonal algorithm with second-order implicit smoothing, and (3) the diagonal algorithm 
with fourth-order implicit smoothing. The gener 2 d convergence histories are the same for all three 
cases, and they are all able to converge about three-and-a-half orders of magnitude in less than 
1000 iterations. 


Circular Cylinder 

The flow over a two-dimensional circular cylinder is a simple example of external flow 
over a bluff body which will have separation. Some of the flow quantites from the calculation 
were compared to experimental values in reference 1 and agreement was found. 

An 0-grid which wraps around the cylinder is used for the grid. The (- and ?;-directions 
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Figure 1.- Separation length versus Reynolds number for flow over a backward facing step. 



ITERATION NUMBER 

Figure 2.- Convergence history for flow over a backward facing step. 
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are shown in Figure 3. The ( axis points in the radial direction and the rj axis points along the 
circumferential direction since any periodic direction must be associated with the rj direction. 


FLOW 



Figure 3.- Schematic diagram of an 0-grid around a circular cylinder. 

The relevant parameters used in this calculation are given in the following list. 


BETA 

= 5 

KPERI 

= 1 

DTAU 

= 0.1 

LMAX 

= 41 

ENDACC 

= 1 

SMU 

= 0.1 

lORTHO 

= 2 

SMUIM 

= 0.3 

JMAX 

KMAX 

= 3 
= 80 

SMUPRS 

= 0.1 


The convergence history is shown in Figure 4 in which the log of RMSDQ and the log of 
RMSDIV are plotted versus the iteration number. The variable RMSDQ is the rms of the change 
of all the variables during one iteration and RMSDIV is the rms of the divergence of velocity 
calculated during that time-step. 

This calculation is typical of most external flow problems in that the calculation is not as 
sensitive to the downstream boundary conditions as in an intern 2 d flow problem. This calculation 
was found to work well when the downstream variables were updated using only a second-order 
extrapolation. 


Three-Dimensional Post 


The geometry of this calculation consists of a circular post mounted between two parallel 
flat plates, such that the post is perpendicular to the plates and the incoming flow is parallel to 
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the two plates. The calculations used one of the flat plates in the xy plane at z = 0 and a 
symmetry boundary condition in the xy plane at z = 5 post diameters. For this flow calculation, 
a series of 29 identical 35-by-41 0-grids similar to the one described for the two-dimensional 
circular cylinder were used parallel to the flat plate, with clustering close to the flat plate. The 
inflow was prescribed to be a partially-developed boundary layer flow between the flat plates. 

The relevant parameters used in this calculation are given in the following list. 


BETA 

= 5 

KPERI 

= 1 

DTAU 

= 0.05 

LMAX 

= 29 

ENDACC 

= 1 

SMU 

= 0.1 

lORTHO 

= 2 

SMUIM 

= 0.3 

JMAX 

KMAX 

= 35 
= 41 

SMUPRS 

= 0.1 


In Figure 5, the geometry is shown along with the velocity vectors denoting the inflow 
velocity profile. The periodic boundary can be seen on the downstream side of the post. 

The program output for this example is included here. The first item printed to the 
output includes the veilues of all the variables in namelist DATAIN. Next, the convergence history 
is printed every five iterations. Note that the format of the output was changed here from that in 
the code in order to fit it on the page. The entry in the first column is NT, which is the current 
iteration number. The second entry, RMSDQ, is the rms of the change in all four variables at all 
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Figure 5.- Three-dimensional post mounted normal to a flat plate with parabolic inlet velocity 
profile. 

the grid points. The entry RMSCO is the rms of the right-hand side of the continuity equation. 
This includes the divergence of the velocity field plus the smoothing terms. The fourth entry 
is RMSDIV, which is the rms of the divergence of the velocity field. The entry DQMAX is the 
maximum change of any of the variables occurring in that iteration, and the J, K, and L values 
give the location of the maximum change. 

Both RMSDQ and DQMAX are measures of convergence. A solution is considered 
converged if these quantities decrease by three to four orders of magnitude over their original 
values. In this run, RMSDQ and DQMAX have both converged over three orders in 600 iterations. 
RMSDIV is a measure of how well the velocity field satisfies the divergence-free condition. This 
code can usually reduce RMSDIV to between 10“^ and 10~®. To achieve this reduction in the 
present case, it was necessary to reduce the amount of explicit smoothing added to the continuity 
equation. To do this, one can just change the value of SMUARY(l) in subroutine SMOOTH. This 
is initially set equal to the input variable SMUPRS. In this calculation, SMUPRS is initially set 
equal to 0.5. It can be seen that as the calculations progress, RMSCO continues to converge while 
RMSDIV approaches a constant. This indicates that the smoothing is dominating the continuity 
equation and should be reduced. The value of SMUARY(l) was decreased by a factor of two at 
time steps of 350 and 450. The net effect was to decrease the divergence further. It should be 
noted that if the smoothing is decreased too much, the pressure field will start to develop kinks 
and the calculation could become unstable. 
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&DATAIN DISKOUT = .T., PLOT3D = .F., JMAX = 35, KMAX = 41, LMAX = 29, NTMAX 
= 600, NPRNT = 600, IPRNT = 5, TIMACC = 1, ENDACC = 1, DTAU = 5.E-2, BETA = 5., 
SMU = 0.1, SMUIM = 0.3, SMUPRS = 0.5, REYNUM = 100., DXDT = 0., RLl = 1., KPERI 
= 1, ISTART = 0, IBLKDIA = 2, IMPSMO = 2, &END 


NT 

RMSDQ 

RMSCO 

RMSDIV 

DQMAX 

J 

K 

L 

5 

0.2590E-01 

0.6427E-01 

0.1865E+00 

0.3160E+00 

2 

22 

21 

10 

0.2643E-01 

0.7407E-01 

0.2856E+00 

-0.3742E+00 

2 

12 

19 

15 

0.2471E-01 

0.6486E-01 

0.2394E+00 

0.3503E+00 

2 

35 

15 

20 

0.2211E-01 

0.5559E-01 

0.2212E+00 

0.4163E+00 

2 

11 

16 

25 

0.1747E-01 

0.4495E-01 

0.1856E+00 

0.2632E+00 

2 

30 

18 

30 

0.1359E-01 

0.3389E-01 

0.1383E+00 

0.1589E+00 

2 

29 

28 

35 

0.1052E-01 

0.2561E-01 

0.9981E-01 

-0.1227E+00 

2 

13 

16 

40 

0.7825E-02 

0.1923E-01 

0.7474E-01 

0.8304E-01 

6 

21 

20 

45 

0.5990E-02 

0.1467E-01 

0.6023E-01 

0.6616E-01 

2 

31 

15 

50 

0.4737E-02 

0.1146E-01 

0.5052E-01 

0.5518E-01 

2 

14 

17 

55 

0.3679E-02 

0.8812E-02 

0.3993E-01 

0.4342E-01 

10 

22 

28 

60 

0.2918E-02 

0.6932E-02 

0.3002E-01 

0.3717E-01 

11 

22 

28 

65 

0.2343E-02 

0.5345E-02 

0.2278E-01 

0.3016E-01 

12 

22 

28 

70 

0.1897E-02 

0.4236E-02 

0.1959E-01 

0.2410E-01 

13 

22 

28 

75 

0.1594E-02 

0.3451E-02 

0.1788E-01 

0.1939E-01 

14 

22 

28 

80 

0.1333E-02 

0.2733E-02 

0.1578E-01 

0.1584E-01 

15 

22 

28 

85 

0.1121E-02 

0.2210E-02 

0.1336E-01 

0.1310E-01 

16 

21 

28 

90 

0.9644E-03 

0.1783E-02 

0.1145E-01 

0.1083E-01 

17 

21 

28 

95 

0.8348E-03 

0.1445E-02 

0.1075E-01 

0.8746E-02 

18 

21 

28 

100 

0.7383E-03 

0.1228E-02 

0.1066E-01 

-0.6955E-02 

7 

38 

28 

105 

0.6612E-03 

0.1057E-02 

0.1036E-01 

-0.6311E-02 

7 

38 

28 

110 

0.5952E-03 

0.9153E-03 

0.9791E-02 

-0.5887E-02 

8 

38 

28 

115 

0.5403E-03 

0.8070E-03 

0.9315E-02 

-0.5473E-02 

8 

38 

28 

120 

0.4921E-03 

0.7120E-03 

0.9122E-02 

-0.5041E-02 

8 

38 

28 

125 

0.4506E-03 

0.6385E-03 

0.9082E-02 

-0.4583E-02 

8 

38 

28 

130 

0.4153E-03 

0.5804E-03 

0.8995E-02 

-0.4113E-02 

8 

38 

28 

135 

0.3844E-03 

0.5278E-03 

0.8830E-02 

-0.3678E-02 

9 

38 

28 

140 

0.3570E-03 

0.4825E-03 

0.8675E-02 

-0.3334E-02 

9 

38 

28 

145 

0.3326E-03 

0.4416E-03 

0.8581E-02 

-0.3005E-02 

9 

38 

28 

150 

0.3114E-03 

0.4064E-03 

0.8525E-02 

-0.2773E-02 

16 

1 

28 

155 

0.2933E-03 

0.3792E-03 

0.8466E-02 

-0.2560E-02 

16 

1 

28 

160 

0.2784E-03 

0.3594E-03 

0.8394E-02 

-0.2364E-02 

17 

1 

28 

165 

0.2659E-03 

0.3459E-03 

0.8329E-02 

-0.2197E-02 

17 

1 

28 

170 

0.2554E-03 

0.3369E-03 

0.8280E-02 

-0.2032E-02 

17 

1 

28 

175 

0.2466E-03 

0.3321E-03 

0.8241E-02 

-0.1868E-02 

17 

1 

28 

180 

0.2394E-03 

0.3318E-03 

0.8202E-02 

-0.1722E-02 

18 

1 

28 

185 

0.2339E-03 

0.3357E-03 

0.8162E-02 

-0.1584E-02 

18 

1 

28 

190 

0.2302E-03 

0.3435E-03 

0.8125E-02 

-0.1449E-02 

18 

1 

28 

195 

0.2283E-03 

0.3546E-03 

0.8094E-02 

-0.1330E-02 

19 

1 

28 

200 

0.2282E-03 

0.3686E-03 

0.8068E-02 

-0.1258E-02 

22 

1 

28 

205 

0.2295E-03 

0.3848E-03 

0.8044E-02 

-0.1246E-02 

22 

1 

28 

210 

0.2319E-03 

0.4019E-03 

0.8024E-02 

-0.1201E-02 

21 

1 

28 

215 

0.2345E-03 

0.4176E-03 

0.8008E-02 

-0.1130E-02 

21 

1 

28 

220 

0.2363E-03 

0.4301E-03 

0.7997E-02 

-0.1035E-02 

21 

1 

28 

225 

0.2366E-03 

0.4375E-03 

0.7990E-02 

-0.9277E-03 

20 

1 

28 

230 . 

0.2349E-03 

0.4392E-03 

0.7984E-02 

-0.8438E-03 

2 

8 

27 

235 

0.2312E-03 

0.4352E-03 

0.7982E-02 

-0.8612E-03 

2 

8 

27 
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NT 

RMSDQ 

RMS CO 

RMSDIV 

DQMAX 

J 

K 

L 

240 

0.2259E-03 

0.4266E-03 

0.7983E-02 

-0.8738E-03 

2 

8 

27 

245 

0.2194E-03 

0.4148E-03 

0.7989E-02 

-0.8813E-03 

2 

8 

28 

250 

0.2124E-03 

0.4014E-03 

0.7999E-02 

-0.8814E-03 

2 

8 

28 

255 

0.2054E-03 

0.3875E-03 

0.8014E-02 

-0.8692E-03 

2 

8 

28 

260 

0.1987E-03 

0.3742E-03 

0.8033E-02 

-0.8422E-03 

2 

8 

28 

265 

0.1926E-03 

0.3621E-03 

0.8056E-02 

-0.8018E-03 

3 

9 

28 

270 

0.1871E-03 

0.3513E-03 

0.8081E-02 

-0.7584E-03 

3 

34 

28 

275 

0.1821E-03 

0.3416E-03 

0.8107E-02 

-0.7067E-03 

3 

34 

28 

280 

0.1774E-03 

0.3324E-03 

0.8133E-02 

-0.6500E-03 

3 

34 

28 

285 

0.1730E-03 

0.3235E-03 

0.8157E-02 

-0.6020E-03 

6 

39 

28 

290 

0.1685E-03 

0.3143E-03 

0.8180E-02 

-0.5649E-03 

6 

39 

28 

295 

0.1639E-03 

0.3047E-03 

0.8200E-02 

-0.5304E-03 

7 

39 

28 

300 

0.1590E-03 

0.2943E-03 

0.8218E-02 

-0.5028E-03 

7 

39 

28 

305 

0.1539E-03 

0.2833E-03 

0.8233E-02 

-0.4749E-03 

7 

39 

28 

310 

0.1484E-03 

0.2717E-03 

0.8246E-02 

-0.4461E-03 

7 

39 

28 

315 

0.1428E-03 

0.2595E-03 

0.8257E-02 

-0.4187E-03 

8 

39 

28 

320 

0.1371E-03 

0.2470E-03 

0.8265E-02 

-0.4151E-03 

34 

1 

28 

325 

0.1314E-03 

0.2345E-03 

0.8273E-02 

-0.4132E-03 

34 

1 

28 

330 

0.1259E-03 

0.2223E-03 

0.8279E-02 

-0.4127E-03 

34 

1 

28 

335 

0.1206E-03 

0.2105E-03 

0.8285E-02 

-0.4134E-03 

34 

1 

28 

340 

0.1156E-03 

0.1993E-03 

0.8291E-02 

-0.4148E-03 

34 

1 

28 

345 

0.1108E-03 

0.1888E-03 

0.8296E-02 

-0.4160E-03 

34 

1 

28 

350 

0.2740E-03 

0.1046E-02 

0.8301E-02 

0.5572E-02 

2 

27 

24 

355 

0.1444E-03 

0.3286E-03 

0.4654E-02 

-0.1646E-02 

3 

13 

19 

360 

0.1197E-03 

0.2430E-03 

0.4369E-02 

0.9543E-03 

34 

2i 

4 

365 

0.1081E-03 

0.2149E-03 

0.4534E-02 

-0.1044E-02 

2 

21 

15 

370 

0.1031E-03 

0.1923E-03 

0.4494E-02 

-0.6718E-03 

2 

21 

28 

375 

0.9620E-04 

0.1697E-03 

0.4496E-02 

-0.5160E-03 

7 

21 

19 

380 

0.8872E-04 

0.1475E-03 

0.4491E-02 

-0.4101E-03 

9 

21 

20 

385 

0.8326E-04 

0.1302E-03 

0.4464E-02 

-0.3655E-03 

34 

1 

28 

390 

0.7902E-04 

0.1187E-03 

0.4466E-02 

-0.3348E-03 

34 

1 

28 

395 

0.7593E-04 

0.1091E-03 

0.4488E-02 

-0.3161E-03 

34 

1 

2 

400 

0.7371E-04 

0.1037E-03 

0.4496E-02 

-0.3094E-03 

34 

1 

3 

405 

0.7217E-04 

O.lOllE-03 

0.4487E-02 

0.2999E-03 

34 

41 

23 

410 

0.7105E-04 

0.9960E-04 

0.4480E-02 

0.3009E-03 

34 

41 

23 

415 

0.7020E-04 

0.9798E-04 

0.4485E-02 

0.3018E-03 

34 

1 

23 

420 

0.6984E-04 

0.9651E-04 

0.4496E-02 

0.3027E-03 

34 

1 

23 

425 

0.7039E-04 

0.9653E-04 

0.4507E-02 

0.3030E-03 

34 

1 

23 

430 

0.7201E-04 

0.9873E-04 

0.4516E-02 

0.3025E-03 

34 

1 

23 

435 

0.7402E-04 

0.1018E-03 

0.4524E-02 

0.3230E-03 

2 

22 

19 

440 

0.7520E-04 

0.1029E-03 

0.4532E-02 

0.3490E-03 

2 

21 

28 

445 

0.7460E-04 

O.lOOlE-03 

0.4536E-02 

0.3544E-03 

4 

34 

17 

450 

0.1574E-03 

0.5739E-03 

0.4537E-02 

0.3060E-02 

3 

27 

23 

455 

0.9444E-04 

0.2107E-03 

0.2630E-02 

0.1044E-02 

2 

22 

23 

460 

0.7930E-04 

0.1491E-03 

0.2430E-02 

0.6952E-03 

4 

21 

28 

465 

0.6833E-04 

0.1174E-03 

0.2499E-02 

-0.5334E-03 

2 

21 

16 

470 

0.6339E-04 

0.1012E-03 

0.2459E-02 

-0.3949E-03 

2 

21 

27 

475 

0.6121E-04 

0.9602E-04 

0.2449E-02 

-0.3006E-03 

7 

21 

20 

480 

0.5949E-04 

0.9406E-04 

0.2451E-02 

0.2932E-03 

34 

1 

24 

485 

0.5911E-04 

0.9624E-04 

0.2441E-02 

0.2958E-03 

34 

1 

24 

490 

0.5994E-04 

0.1017E-03 

0.2441E-02 

0.2991E-03 

34 

1 

24 
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NT 

RMSDQ 

RMSCO 

RMSDIV 

DQMAX 

J 

K 

L 

495 

0.6120E-04 

0.1064E-03 

0.2452E-02 

0.3032E-03 

34 

1 

25 

500 

0.6211E-04 

0.1105E-03 

0.2457E-02 

0.3084E-03 

33 

1 

26 

505 

0.6255E-04 

0.1135E-03 

0.2453E-02 

0.3124E-03 

33 

1 

26 

510 

0.6240E-04 

0.1150E-03 

0.2448E-02 

0.3146E-03 

34 

1 

26 

515 

0.6119E-04 

0.1136E-03 

0.2448E-02 

0.3161E-03 

34 

1 

26 

520 

0.5852E-04 

0.1080E-03 

0.2449E-02 

0.3156E-03 

34 

1 

28 

525 

0.5452E-04 

0.9839E-04 

0.2449E-02 

0.3124E-03 

34 

1 

28 

530 

0.4991E-04 

0.8596E-04 

0.2448E-02 

0.3062E-03 

34 

1 

28 

535 

0.4572E-04 

0.7300E-04 

0.2449E-02 

0.2968E-03 

34 

1 

28 

540 

0.4268E-04 

0.6205E-04 

0.2452E-02 

0.2848E-03 

34 

1 

28 

545 

0.4083E-04 

0.5468E-04 

0.2455E-02 

0.2710E-03 

34 

1 

28 

550 

0.3970E-04 

0.5076E-04 

0.2458E-02 

0.2563E-03 

34 

1 

28 

555 

0.3876E-04 

0.4909E-04 

0.2459E-02 

0.2417E-03 

34 

1 

28 

560 

0.3782E-04 

0.4871E-04 

0.2460E-02 

0.2282E-03 

34 

1 

28 

565 

0.3696E-04 

0.4933E-04 

0.2460E-02 

0.2162E-03 

34 

1 

28 

570 

0.3630E-04 

0.5075E-04 

0.2461E-02 

0.2064E-03 

34 

1 

28 

575 

0.3585E-04 

0.5254E-04 

0.2462E-02 

0.1987E-03 

34 

1 

28 

580 

0.3550E-04 

0.5415E-04 

0.2463E-02 

0.1933E-03 

34 

1 

28 

585 

0.3512E-04 

0.5515E-04 

0.2464E-02 

0.1899E-03 

34 

1 

28 

590 

0.3461E-04 

0.5531E-04 

0.2464E-02 

0.1884E-03 

34 

1 

28 

595 

0.3390E-04 

0.5461E-04 

0.2465E-02 

0.1885E-03 

34 

1 

28 

600 

0.3299E-04 

0.5311E-04 

0.2465E-02 

0.1898E-03 

34 

1 

28 


20 



REFERENCES 


1. Kwak, D.; Chang, J. L. C.; Shanks, S. P.; and Chakravarthy, S.: A Three-Dimensional 
Incompressible Navier-Stokes Flow Solver Using Primitive Variables. AIAA J., vol 24, 
no. 3, 390-396, Mar. 1986 

2. Chorin, A. J.: A Numerical Method for Solving Incompressible Viscous Flow Problems. 
J. Comp. Phys., vol. 2, 1967, pp. 12-26. 

3. Steger, J. L.; and Kutler, P.: Implicit Finite- Difference Procedures for the Computation 
of Vortex Wakes. AIAA J., vol. 15, no. 4, 1977, pp. 581-590. 

4. Chang, J. L. C.; and Kwak, D.: On the Method of Pseudo Compressibility for Numeri- 
cally Solving Incompressible Flows. AIAA Paper 84-0252, 1984. 

5. Beam, R. M.; and Warming, R. F.: An Implicit Finite-Difference Algorithm for Hyper- 
bolic Systems in Conservation-Law Form. J. Comp. Phys., vol. 22, 1976, pp. 87-110. 

6. Briley, W. R.; and McDonald, H.: Solution of the Three-dimensional Compressible 
Navier-Stokes Equations by an Implicit Technique. Proceedings of the Fourth Interna- 
tional Conference on Numerical Methods in Fluid Dynamics, Lecture Notes in Physics, 
vol. 35, Springer- Verlag, New York 1975, pp. 105-110. 

7. Rogers, S. E.; Chang, J. L. C.; and Kwak, D: A Diagonal Algorithm for the Method of 
Pseudocompressibility. AIAA Paper 86-1060, 1986. 

8. Rogers, S. E.; Kwak, D.; and Kaul, U.: On the Accuracy of the Pseudocompressibility 
Method in Solving the Incompressible Navier-Stokes Equations. Appl. Math. Modelling, 
vol. 11, 1987, pp. 35-44. 

9. Zilliac, G.: A Computational/Experimental Study of the Flow Around a Body of Rev- 
olution at Angle of Attack. NASA TM-88329, 1986. 

10. Champney, J. M.; Coakley, T. J.; and Kaul, U. K.: Incompressible Viscous Flow Simu- 
lations of the NFAC Wind Tunnel. ATM-FR-25-012, 1986. 

11. Kaul, U. K.; Kwak, D.; and Wagner, C: A Computational Study of Saddle Point Sepa- 
ration and Horseshoe Vortex System. AIAA Paper 85-0182, 1985. 

12. Rogers, S. E.; Kwak, D.; and Kaul, U. K.: A Numerical Study of Three-Dimensionsal 
Incompressible Flow Around Multiple Posts. AIAA Paper 86-0353, 1986. 

13. Chang, J. L. C.; Kwak, D.; and Dao, S. C.: A Three Dimensional Incompressible Flow 
Simulation Method and its Application to the Space Shuttle Main Engine, Part 1, Lam- 
inar Flow. AIAA Paper 85-0175, 1985. 


21 



14. Yang, R. J.; Chang, J. L. C.; and Kwak, D.: A Navier-Stokes Simulation of the Space 
Shuttle Main Engin Hot Gas Manifold. AIAA Paper 87-0368, 1987. 

15. Armaly, B. F.; Durst, F.; Pereira, J. C. F.; and Schonung, B.: Experimented and Theo- 
retical Investigation of Backward Facing Step Flow. J. Fluid Mech., vol. 127, 1983, pp. 
473-496. 

16. Warming, R. F.; Beam, R. M.; and Hyett, B. J.: Diagonalization and Simultaneous 
Symmetrization of the Gas-Dynamic Matrices. Mathematics of Computation, vol. 29, 
1975, pp. 1037-1045. 

17. Turkel, E.: Symmetrization of the Fluid Dynamic Matrices with Applications, Mathe- 
matics of Computation, vol. 27, 1973, pp. 729-736. 

18. Pulliam, T. H.; and Chaussee, D. S.: A Diagonal Form of an Implicit Approximate- 
Factorization Algorithm. J. Comp. Phys., vol. 39, 1981, pp. 347-363. 

19. Flores, J.: Convergence Acceleration for a Three-Dimensional Euler/Navier-Stokes Zoned 
Approach. AIAA Paper 85-1495, 1985. 


22 



APPENDIX 

THEORETICAL FOUNDATION 


I. Governing Equations 

Incompressible Navier-Stokes Equations 

Unsteady, three-dimensional, viscous, incompressible flow with constant density is gov- 
erned by the following Navier-Stokes equations: 


dui 


+ 


dui 

dxi 

dutUj 


= 0 


dp , dr, 


+ 


( 1 ) 


( 2 ) 


dt dxj dxi dxj 

where t is the time, Xi the Cartesian coordinates, the corresponding velocity components, p 
the pressure, and Tij the viscous stress tensor. All of the variables are nondimensionalized by a 
reference velocity and length scale. The viscous stress tensor can be written as 


T{j — luSij Rij 


( 3 ) 

Here, u is the kinematic viscosity, Sij is the strain rate tensor, and Rij is the Reynolds 
stresses. The strain rate tensor is defined by 


\ dui duj 
2^dxj^ dxi^ 


( 4 ) 


Various levels of closure models for Rij are possible. In the present code, turbulence is 
simulated by an eddy viscosity model, using a constitutive equation of the following form: 

Rij = -Rkk^ij — 2utSij (5) 

where ut is the turbulent eddy viscosity. By including the normal stress, Rkk, in the pressure, 
the variable u in equation (3) can be replaced by (u -|- ut). 


Tij = 2(v + ut)Sij = 2urSij ( 6 ) 

In the remainder of this paper, the totad viscosity, ut^ will be represented simply by u. Therefore, 
the incompressible Navier-Stokes equations are modified to allow variable viscosity in the present 
formulation. 


Pseudocompressible Formulation 

The system of equations given in the preceding section is not as readily solvable by a fast 
alternating-direction-implicit (ADI) scheme as the unsteady compressible Navier-Stokes equations 
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which are hyperbolic. This difficulty can be alleviated by modifying the continuity equation as 
follows: 


^ ^ 0 
(3 dt dxi 


( 7 ) 


By combining equations (2) and (7), the governing equations are written as one vector 


equation 


where 


|d+A(£-E.) + |(F-F.)+A(g-G.) = 0 


D = 


p 

u 

V 

LtUJ 


( 8 ) 



1 — 
3 


■ ■ 





vr 

F = 

vu 

(? = 

wu 


uv 


+p 


wv 


- uw . 


. VV) . 


kj 

+ 

1 


( 9 ) 


1 

• 0 ■ 

1 1 

• 0 ■ 


- 0 ■ 

II 

'^xy 

II 

Ty* 

Tyy 


'^ZX 

1 

- '^xz - 

1 1 

LryJ 


.Tzz. 


Coordinate Transformation 

To perform calculations on three-dimensioned arbitrarily shaped geometries, the follow- 
ing generalized independent variables are introduced which transform the physical coordinates 
into general curvilinear coordinates. 


T = t 

V = 

C ~ C(®>2/> ^>0 


( 10 ) 


The spatial and time derivatives become 

dxi dxi d( ^ dx{ drj ^ dx{ dC dxi d^j 
d _ d dii d 

di~ dr dt dij 
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Here the Jacobian of the transformation is defined as 


J = 


d(x,y , ») 


^as iz 

Vx Vy Vz 

Cx Cv Cx 


where, for example, 




— aif Vv - a'. 


drj 

dy 


where 


In actual coding, metric terms are calculated as follows : 

1 / VvH - y<^v 

— I Z^Zr, - Xr,Z( 

y<H - ytH 
Hyi - ^iy< 


yi^r, - vnH 

= T, I 



^iy-n - 


J'= det 




Xr, X( 

yi yv y< 


Applying the transformation to equation (8) yields 

_ £.) + - F.) + ^(G - G.) = 0 


where 



E = j[^tD + (,E + iyF + (,G] 
F = — \rjtD + rj^E + rfyF + »/*G] 

G = i[<,Z) + ai; + C,f' + CG] 

Ev = -j\ ixEv + iyFv + 

^t. = J [ rjxEv + riyFy + ’n^G^] 

Gv = j[ CxE,+CyF,+CG,] 


( 12 ) 


( 13 ) 


( 14 ) 


( 15 ) 
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The contravariant velocities, U., V, and W without metric normalization are defined as 


U = (t + + ^yV + izW 

V = »/t + + r}tW 

^ = Ct + Cxtt + ^yV + 


The flux vectors can be rewritten as 


E = 


F = 


G = 


1 

J 


1 

J 


1 

J 


r^U + Up-fiY 

' uU^-Up 

VU + (yP 
tvU + (tp 

^I3V + rjtip-^y 
uV + r}t,p 
vV + 'HyP 

+ rizp 

r^W + Up-P) 

uW + CxP 
vW + CyP 
wW + CzP 


( 16 ) 


(IT) 


Viscous Terms 


The viscous terms are quite lengthy and therefore are given here separately. 

d 


dxj 


Tij = 


9 „ 


j9 

dx 


where, for example, 


«x +«x 
V* +«y 
tWx + Uz 


u. 


+ 


dy 


Uy + V* 

Vy + Vy 
Wy + Vj 


d 


Uz + Wz 
Vz + Wy 
Wz+Wz 


du 

di' 


Uy = 


du 

dy 


When V is constant, the contribution of the second terms in each bracket sums up to zero when 
the velocity field is divergence-free. However, since in general, v varies in space and time these 
terms have to be kept. Then the viscous terms in transformed coordinates are given by 





-0- 

U 

X du dv dw^ 

1 


V 

+ Ux^ + ^V^+^,^) 

'By^' 





1 

•««» 

1 

> 
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A ^ I 

Fv--j{ 


(V^-(V^^ + Vr;- + VC^) 


Ol 

u 

V 

Ltl>J 


. du dv dw ^ 


A- 

f!: 

m^i\ 




ro 

u 

V 

WJ 


,. dv, . dv . dw^ 


A. 


} (18a) 


The viscous terms given in equation (18a) can be simplified under special conditions. If an 
orthogonal coordinate system is used, then 


= 0 for t ^ j 


Equation (18a) becomes 




dU ^*^dCi ^‘dU' 




d(i '^dii '^dCi' 


G. = ^\ (a^ +Cv^ + + (C.~ +c,~ 




dw 


0 

9 f. 

d t. 

dy^* 

o>V 

'dx^i 

L d^U J J 

d^^i 

L J J 


where 


Jm = 


/ 


0 0 0 01 
0 10 0 
0 0 10 
LO 0 0 IJ 


(186) 


For constant v (i.e., laminar flow), the contribution by the second term in parentheses in equation 
(18b) will approach zero as the steady state is reached, where equation (1) is numerically satisfied. 
Therefore, for flow with constant v in orthogonal coordinates, equation (18b) can be reduced 
further as follows: 
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(18c) 


= (7) = 72^ 

Gw = (y) (c®^ + Cy^ + — 73-D 


II. Numerical Algorithm 


Time- Advancing Algorithm 

The numerical algorithm used to advance equation (15) in time is an implicit, nonit- 
erative, approximately factored, finite-difference scheme described by both Beam and Warming 
(ref. 5) and Briley and McDonald (ref. 6). The time-differencing used by this scheme, generally 
known as the trapezoidal rule, is given by 


7,n+i = pn 


Ar 






-k O(Ar’) 


(19) 


where the superscript n refers to the n*^ time-step. The finite-difference form of equation (15) is 

Srb + S^{E - E^) + «,(F - Fw) + ^c(G - Gw) = 0 (20) 

where, for example, 8^ is the finite difference operator for Expressions for F, F,and G 

are given by equation (17); expresions for Fw, Fw,and Gv are given by equation (18a) for a 
nonorthogonal grid, by equation (18b) when an orthogonal grid is assumed, or by equation (18c) 
for the case of laminar flow and an orthogonal grid. By substituting equation (20) into equation 
(19) and using D — DJ^ one obtains 


[SdE - Fw)""^^ + Sr,(F - Fw)’^+^ + 8({G - Gw)”+*] 

Tn-f 1 A'r ^ ^ A ^ 

= J- J’'+‘ [<((£ - - F.r + S((G - G.)“] 


(21) 


Solving for is nonlinear in nature since F”"*"* = E(D'*'*‘^) is a nonlinear function 

of as are F”"*"^ and The following linearization procedure is used. A local Taylor 

expansion about yields 


^n+i = ^ i«(F”+* - F”) + O(At^) 

pn+l ^npn-M _ ^ 0 (At*) (22) 

^n-hl = (jn ^ C'«(F”+* - F”) -h G(At^) 


28 


where Ay and C are the Jacobian matrices 


A 


dE 

dD^ 


. dF 

^ dD'^~dD 


(23) 


The Jacobian matrices can all be represented by the following, where Ai — Ay By or C for i = 
1,2, or 3, respectively. 


Lq -^ 2/3 

L\ Q + Xjti L^v, L%u 

L 2 L\V Q L 2 V X3V 
L 3 L\w L 2 W Q + L 3 W. 

where 

Q — Lq + Liu + L 2 V + X3W 

lo = ({i)«, £i = ((i)., U = (fi)„ £j = («i). 

Vi or ^ for Ai = Ay By or Cy respectively 



(24) 


Substituting equation (22) into equation (21) results in the governing equation in delta form 


|j + [s^iA'^ - Ti) + 6y,{B^ - Tj) + S({C^ - T3)] I (!>’»+' - D^) 

= - Ar |^^(E - Ey,)^ + Sr,{F - Fy,)^ + 6({G - (?v)~j (25) 



where 


riD'‘+^ = = F^+'y r3i>"+* = (?”+* 


h = At for trapezoidal scheme 


= 2At for Euler scheme 


At this point it should be noted that the notation of the form [^((A — r)]D refers to -^(AD) — 
■§^{TD) and not - ^D. 


Approximate Factorization 

The solution of equation (25) would involve a formidable matrix-inversion problem. With 
the use of an ADI-type scheme, the problem could be reduced to the inversion of three matrices 
of smeill bandwidth, for which there exists efficient solution algorithms. The particular ADI form 
used here is known as approximate factorization (AF) (refs. 5,6). It is difficult to implement 
the AF scheme to solve equation (25) in its full matrix form. Noting that at steady-state, the 
left-hand side of equation (25) approaches zero, a simplified expression for the viscous term as 
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I 

I 

shown in equation (18c) is used on the left-hand side. To maintain the accuracy of the solution, 
the entire viscous term is used on the right-hand side. Using this, the governing equation becomes 

- D'^) = RHS (26) 

where 

£,= (27) 

and RHS is the right-hand side of equation (25). When second-order central differencing is 
used, the solution to this problem becomes the inversion of three block-tridiagonal matrices. The 
inversion problem is reduced to the three inversions 

(X,)A5 = RHS 

(Z;^)A5 = A.D (28) 

j (X^)AI>"+^ = AD. 

These inversions are carried out for all interior points, and the boundary conditions 
can be implemented explicitly. However, it is possible to implement the boundary conditions 
implicitly, which will be discussed in a later section. 

The factorization has introduced the following second-order cross-product term into the 

equation. 

[S(ASr,B + Sr,BS(C -b 6(C6(A] AD + 0{h^) (29) 

where 

A = A^-^u B = B’‘-72, C7 = C^-73, h = ^ 

I To maintain the second-order accuracy of the scheme, the added terms must be kept smaller 

than the original terms in the equations everywhere in the computational domain. This puts a 
restriction on the size of the artificial compressibility parameter /3. The proper choice of yS will be 
; discussed in a later section. When the approximate factorization scheme is applied, the stability 

of the scheme is dependent on the use of higher-order smoothing terms. These are used to damp 
out higher frequency oscillations which arise in the solution because of the use of second-order 
central differencing, and will be discussed in a following section. 

j Numerical Dissipation / Smoothing 

Higher-order smoothing terms are required to make the present algorithm stable. These 
terms help to damp out the higher order oscillations in the solution that are caused by the use of 
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second-order centred differencing. The smoothing term is derived in conjunction with an upwind 
finite-difference approximation. The fourth-order smoothing operator is given by 


(VxA*)^tt = 1 -f 6uj - 4uj+i -f Uj+ 2 ) 

The second-order smoothing operator is given by 

VxAx« = 2 ^^^^+^ ~ 

By including these smoothing terms, equations (26) and (27) become 


(30) 

(31) 


L(L,L((D'+' - D") = RHS - €.[(V<A«)“ + (V,A,)" + (V^A^)’]©” ( 26 ') 

where 

Lr,= J-|-^J'‘+'M5"-72) + CiV,A, (27') 

Xc= /+^J”+‘^a<^'‘- 73 ) + eiVcA^ 

Here is the implicit smoothing coefficient and €« is the explicit smoothing coefficient. 

To preserve the tridiagonal nature of the system, only second-order smoothing can be 
used on the implicit side of the equation, whereas fourth-order smoothing is used on the right- 
hand side. However, when the diagonal algorithm is used, the bandwidth of the system can 
be increased to a pentadiagonal, since the scaler inversions are relatively cheap. This makes it 
possible to use fourth-order smoothing on the implicit side of the equations as well. 

When using the second-order implicit, fourth-order explicit scheme, in order to damp 
out the numerical fluctuations as time advances, the smoothing parameter must satisfy 

€i > 2te (32) 

In practice, = 3ce is often used, with c* ^ 0-1 • It sometimes is necessary to change the amount 
of smoothing applied explicitly to the pressure equation; this will be discussed in a later section. 

Diagonal Algorithm 

In a diagonal algorithm, a similarity transform is used to diagonalize the Jacobian ma- 
trices and uncouple the set of equations. The equations can then be solved by solving scalar 
tridiagonal matrices instead of solving block tridiagonal matrices. A similarity transform which 
symmetrizes and diagonalizes the matrices of the compressible gas-dynamic equations has been 
used by Warming et al. (ref. 16) and Turkel (ref. 17). This method was used by Pulliam and 
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Chaussee (ref. 18) to produce a diagonal algorithm for the Euler equations. This can be applied 
to the compressible Navier-Stokes equations to obtain a considerable savings in computing time 
(ref. 19). In this section, similarity transforms for the matrices used in the method of pseudocom- 
pressibility are presented. They are used in a diagonal algorithm which results in a substantial 
reduction in computer time. 

Similarity transformations exist which diagonedize the Jacobian matrices. 

Ai = TikiTr^ (33) 

where Aj is a diagonal matrix whose elements are the eigenvalues of the Jacobian matrices. It is 
given by 


A< = i 


(34) 


■Q 0 0 0 

0 Q 0 0 

0 0 Q+C 0 

.0 0 0 Q-C} 

where C is the pseudospeed of sound, which is given by 

C=^J{Q- LoY + !3{L\ 

The Ti matrix is composed of the eigenvectors of the Jacobian matrix. For i=l (^-sweep), the 
first two eigenvectors are given by 


For i=2 ( 77 -sweep) they are 


For i=3 ((^-sweep) they are 



■ 0 ■ 


■ 0 ■ 

= 

k 

, X2 = 

k 

k 




- k ■ 


• 0 ■ 


■ 0 ■ 

z= 

x^ 

k 

II 

k 

k 




Li<J 


• 0 - 


■ 0 ■ 

= 

k 

II 

X( 

k 


. Zr) . 


J 


(35) 


(36) 


(37) 


The terms in these eigenvectors have been normalized to simplify the inverse matrix with, for 
example, 

- 

“ [L\ + L| + 111*'* 


= 


Xr,VJ 


[Ll + Ll + Llf^ 


32 



The eigenvectors associated with the last two eigenvalues are given by 


C{C-Q)_ 1 
u + Li{C — Q) 
v + Lj{C_-Q) 
.w + Li{C — Q). 


X4 


C{C + Q) 1 
u-Li{C + Q) 
v-L2{C + Q) 
.w-L3{C + Q)\ 


( 38 ) 


where 


Q = Liu + L2V + L3W 

C = ^Q^+fi 

Li 

* y/L\ ^Ll^L\ 


The columns of the Ti matrix are formed with the eigenvectors 


Ti=[Xi X2 Xi X4] 


The elements of the inverse of the eigenvector matrix are given by 

- Z 3 V) + y<.+i(^3U - Liw) 

+ ^(i+x(Liv - L2U)] 

(Tr^ )2i = (LiV - L 2 W) + {Liw - Uu) 

+ %+j(^2« - Liv)] 


{Trhi = 




41 = 


2C72 

1 

2C^ 

1 r- 


{Ti )i 2 = + Qv) - + Qw)] 


(?;"')22 


(^■')32 

irr%2 


+ Qw) - %+,(I-2/3 + Qv)] 


Li(Q + C) 

2C2 

LliQ-C) 

2C^ 
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(Tr'hz 

( r .“^)23 

(Tru 

(Tru 

(T^U 

(T-“^)34 


i [*(<+! (^ 3 /? + Qw) - Z(,^.,(Li/3 + Qu)] 
i[% + 2(^l^ + Q«) - ^+2(^3^ + Q«^)] 

^ 2 (g + c) 

2^2 

X2(Q - C) 

2<?2 

^[ye<+i(-^i^ + Q“) - + Qw)] 

i[«(i+2(-t2/3 +'Qv) - + Q«)l 

Lz{Q + Q) 

2^2 

X3(Q - C) 

2C2 


( 39 ) 


The determinant of Ti is given by 

det{Ti) = 2C’ 

which remains bounded independent of the geometry. For more detail on the derivation of these 
matrices, see ref. 7. 

The implementation of the diagonal scheme involves replacing the Jacobian matrices in 
the implicit operators with the product of the similarity transform matrices and the diagonal 
matrix as given in equation (33). The identity matrix in the implicit operators is replaced by the 
product of the similarity-transform matrix and its inverse. A modification is made to the implicit 
viscous terms by replacing the Im matrix with an identity matrix so that the transformation 
matrices may also be factored out of these terms. This modification implicitly adds an additional 
viscous dissipation term to the pressure. The transformation matrices are now factored out of 
the implicit operators to give 


£< = r {(/+^«(( A (- 7 .)) 27 * 

c-n = r,(/ + ^ J<,(A, - 7 j)|T-> (40) 


where the implicit viscous terms are now given by 

Ti = 

Since the transformation matrices are dependent on the metric quantities, factoring 
them outside the difference operators introduces an error. No modification has been made to 
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the right-hand side of the equation and therefore these linearization errors will only effect the 
convergence path toward the solution, not the steady-state solution. 

The implementation of this algorithm instead of the block algorithm will result in a 
substantial reduction in computational time per iteration because of the decrease in the number 
of operations performed. Additionally, considerably less memory is required to store the elements 
on the left-hand side. This additional memory was used to further vectorize the existing code 
as follows. Since the solution of a tridiagonal block or scalar matrix is recursive, it cannot be 
vectorized for loops which use the current sweep direction as the inner do-loop index. However, if 
a large number of these matrices are passed into the inversion routines at once, then vectorization 
can take place in the ‘nonsweep’ direction. 


Boundary Conditions 

An important part of any computer code is the proper implementation of boundary 
conditions. The code must be capable of handling the several different types of boundaries 
encountered in numerical simulations which include solid-surface, inflow and outflow, and far- 
field boundaries. 

Solid Surface 


At a solid surface boundary, the usu^ no-slip condition is applied. In general, the grid 
point adjacent to the surface will be sufficiently close so that constant pressure normal to the 
surface in the viscous boundary layer can be assumed. For a ^ = constant surface, this can be 
expressed as 



( 41 ) 


This boundary condition can be implemented either explicitly or implicitly. However, the implicit 
implementation will enhance the stability of the code. This can be done during the ^ sweep by 
including the following in the matrix to be inverted. 


where 


JADj,fc,i 4- cADj^k,2 = / 


■-1 

0 

0 

0 - 


'PL=2 — PL=1 ■ 

0 

0 

0 

0 

A 

/ = 

0 

0 

0 

0 

0 

0 

. 0 

0 

0 

0 . 


0 


( 42 ) 


Inflow, Outflow and Far-field Conditions 

The inflow and outflow boundary conditions for an internal flow problem and the far- 
field boundary conditions for an external flow problem can be handled in much the same way. 
The incoming flow for both problems can be set to some appropriate constant as dictated by 
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the problem. For example, at the inlet to a pipe, the pressure can be set to a constant and the 
velocity profile can be specified. The conditions at the downstream however, are the most difficult 
to provide. Simple upwind extrapolation is not well-posed. The best convergence rate is obtained 
if global mass is conserved. So to ensure the best results, the velocities and pressure are first 
updated using a second-order upwind extrapolation. For an exit at L = LMAX, this is 


Q?max — 


Q n+1 f 

lmax — 1 yAzi J V/ma®— 2 

Azi ^ 


(43) 


where 


Azi — Zlmax ^lmax—1 

Az2 = Zimax ^lmax—2 


Then these extrapolated velocities are integrated over the exit boundary to obtain the outlet mass 
flux 

Tnout = I V"' • da. (44) 

J exit 

Then the velocity components are weighted by the mass flux ratio to conserve global mass 


yn+l „ 


^out 


(45) 


If nothing further is done to update the boundary pressure, this procedure can lead to discon- 
tinuities in the pressure because momentum is not being conserved. A method of weighting 
the pressure by a momentum correction was presented by Chang et al (ref. 13). They use the 
following relation to obtain a pressure which corresponds to the mass-weighted velocities 


P 


=p^-~ [(u)VF)’^+^ - (u>W^)*] -h ^(VC ■ VO 


0 


0 


(I)"' - ( 


dc) 


(46) 


where W is the contravariant velocity. In obtaining this formula, it ha^ been assumed that 
the streamlines near the exit plane are nearly straight. Any appreciable deviation will cause a 
discontinuity in the pressure and may lead to an instability. To avoid this, a momentum-weighted 
pressure was used. This was obtained by integrating the momentum-corrected pressure and 
the extrapolated pressure p'^ across the exit 


jn+l 


Jexit 




Ip 


= f p^dd. 

J exit 


The final outlet pressure is then taken to be 


= 


I ‘‘‘P \ 

■f? / 


(47) 


( 48 ) 
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Using these downstream boundary conditions, global conservation of mass and momentum are 
ensured and the scheme will not introduce instabilities into the flow field. 


III. PSEUDOCOMPRESSIBILITY 

In an incompressible flow, a disturbance in the pressure causes waves which travel with 
infinite speed. Introduction of pseudocompressibility causes waves of finite speed, where the 
magnitude of the speed depends upon the pseudocompressibility constant /?. In a true incom- 
pressible flow, the pressure field is affected instantaneously by a disturbance in the flow, but 
with pseudocompressibility, there will be a time lag between the flow disturbance and its effect 
on the pressure field. In viscous flows, the behavior of the boundary layer is very sensitive to 
the streamwise pressure gradient, especially when the boundary layer is separated. If separation 
is present, a pressure wave traveling with finite speed will cause a change in the local pressure 
gradient which will affect the location of the flow separation. This change in separated flow will 
feed back to the pressure field, possibly preventing convergence to a steady state. 


Wave-Vorticity Interactions 


To understand the physical phenomenon of wave propagation and its interaction with 
viscous effect, a simplified one-dimensional form of the governing equation is anlyzed. 


\ dp du - 
du du^ dp 


(49) 


(50) 


where Tv, represents wall shear stress. The pseudospeed of sound for this system of equations is 
given by 

c = 

and the pseudo-Mach number can be defined 


u u 

M = - = - , 

C y/u^ + (3 


(51) 


which is always less than one for any /3 > 0. This equations shows that the present formulation 
will not generate shocks and also that the solution eventually approaches that of incompressible 
flows (see ref. 4 for more detail). In three-dimensional applications, u can be considered to be 
the average velocity in the primary flow direction. 


To investigate the coupling effect between the pressure wave and vorticity, consider a 
channel of length L. The time it takes an upstream propagating wave to travel the distance of 
the computational domain of length L is given by the dimensionless time scale, 

( 52 ) 




c — u 
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The spreading of the viscous effect can be estimated by the boundary-layer growth rate. Therefore, 
for laminar flows, the dimensionless time scale for the viscous effect to spread over a thickness, 
is approximately given by 


Ts 


Re 

4 Xref 


(53) 


To recover the incompressible phenomena, the physics requires that the pressure wave propagate 
much faster than the spreading of the viscous effect, namely. 


Tl « Ts 


(54) 


From this physical reasoning, a general criterion for the lower bound of /3 is obtained given as 
follows: 


/?» 



(55) 


where xs and xx, are the characteristic lengths that the vorticity and the pressure waves have 
to propagate during a given time span. For a channel flow simulation, xx, is equal to the length 
of the channel and xs is hedf the channel width. The computational experiments showed this 
dependency of the pseudocompressibility parameter on the Reynolds number, as well as on the 
characteristic lengths of the geometry (see ref. 8). 


For the near fleld in external flow simulations where pressure waves propagate out to 
infinity, the range of acceptable is less restrictive. The pressure buildup in the near-field region 
is usually temporary, even if (3 is not optimum. However, for internal flows, local pressure buildup 
can be serious when ^ is not properly chosen. ‘ 

The total number of iterations required for the wave to travel downstream and back 
upstream, a total distance of 2T, must be considered. Recalling that the downstream or upstream 
travelling waves propagate with speed (c -f- u) or (c — u) respectively, the total time required for 
one round trip, ti, can be written as 


Tl 


y/u^ +/? 

/3 


2L 


(56a) 


If a constant time step is used throughout the computation, the total number of iterations, iVi, 
for one round trip is 


(566) 

where Ar is the increment of dimensionless time. For convergence, enough computing time, r, 
has to elapse to permit the wave to make at least one round trip. 
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This analysis provides a lower bound for the choice of For the upper bound of the 
error introduced by the approximate factorization will be considered. Recalling equation (29), 
the form of the error is 


E = - r.)«,(5" - r) + . . .] (57) 

This term must be kept smaller than the original terms in the equation. If only the terms which 
contain /3 are included, this restriction can be expressed as 


At 

~T 







or 


At 


2 x; < 1 (58) 

Recalling the expression for given by equation (24), the terms contain (3 give the following 



(59) 


The term to the right of j3 in this inequality is essentially the change in in either the ^,r/,or C 
direction. An estimate of the order of magnitude of this term for the grids generally used is given 
by 



(60) 


which puts the restriction on /? 


0(/3Ar) < 1 


(61) 


For most problems, the restrictions for fi given by equations (56) and (61) are satisfied with a 
value for in the range from 1 to 10. 


Effects of Smoothing on Accuracy 

The explicit smoothing has a large effect on the convergence and accuracy of the pseudo- 
compressibility method. In particular, the explicit smoothing on the pressure can affect whether 
the solution converges to an incompressible flow field. Chang and Kwak (ref. 4) showed that 
the pseudopressure waves decay exponentially with time and vanish as the solution converges; 
thus the change in pressure with time approaches zero. When there is no explicit smoothing 
added to the equation, the divergence of the velocity field identically approaches zero. However, 
when explicit smoothing is included, as the change in pressure approaches zero, the divergence of 
velocity approaches 


dui 

dxi 


[(V<Aj)’ + (V,A,)^ + (V^A^^l p 
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(62) 



where £e(l) is the explicit smoothing parameter for the pressure. If c* is scaled by h, i.e., Cg = 
Ar£e, equation (62) becomes independent of time step. This term can deteriorate the conservation 
of mass depending on the magnitude of /3 and the local pressure gradient. When the pressure 
gradients become substantial, as in the case when a region of separation is present, the smoothing 
term does not approach zero so the divergence of the velocity field is contaminated. In this 
situation, mesh refinement usually does not help in reducing the magnitude of this term. To 
alleviate this problem, it may become necessary to decrease the smoothing coefficient in the 
continuity equation after the solution is almost converged. 
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